function model = projectionStep(vecTheta0,setup)

if setup.dispOn == 1 
    name = 'iter';
else 
    name = 'off';
end
options  = optimset('Display',name,'MaxFunEvals',setup.maxNumEval,...
    'MaxIter',setup.MaxIter,'TolFun',setup.tolF,'TolX',setup.tolX);

if setup.optimizer == 1
    %setupStart = setup;
    %setupStart.xGrid = setupStart.xGrid_2;
    %lambda = 0.0001;
    %tolF    = optimget(options,'TolFun');
    %MaxIterStart = 10;
    %deltaOption = 1;
    %vecThetaOpt = Projection_LMoptimizer(vecTheta0,setupStart,tolF,MaxIterStart,lambda,deltaOption,setup.dispOn);
    
    lambda = 0.10;
    tolF    = optimget(options,'TolFun');
    MaxIter = optimget(options,'MaxIter');
    deltaOption = 0;
    [vecThetaOpt,~,~,exitflag] = Projection_LMoptimizer(vecTheta0,setup,tolF,MaxIter,lambda,deltaOption,setup.dispOn);
elseif setup.optimizer == 2    
    res0 = NLSobjectFunc(vecTheta0,setup);
    Q0   = res0'*res0;
    if isnan(Q0) || isinf(Q0)
        exitflag = -1;
        vecThetaOpt = vecTheta0;
    else
        [vecThetaOpt,~,~,exitflag] = lsqnonlin(@NLSobjectFunc,vecTheta0,[],[],options,setup);
    end
elseif setup.optimizer == 3	
    optionsFmin    = options;
    optionsFmin.MaxIter = 1D4;
    optionsFmin.MaxFunEvals = 1D4;
    [vecThetaOpt,~,exitflag]  = fminsearch(@NLSobjectFunc,vecTheta0,optionsFmin,setup);
elseif setup.optimizer == 4
    [vecThetaOpt,~,exitflag]  = fminunc(@NLSobjectFunc,vecTheta0,options,setup);
elseif setup.optimizer == 5
    % Attempt of an iterative method - does not work well!
    vecTheta = vecTheta0;
    diff = 10;
    maxIter = 2000;
    idx = 1;
    learningRate = 0.05;
    while diff > 1D-4 && idx <= maxIter
        idx = idx + 1;
        controlNew_cu = updateControlsProj(vecTheta,setup);
        [thetaNew,resUnitFree,Yfit] = OLSforProj(controlNew_cu,setup.xGrid,setup.appOrder,setup.scaleX_in_g);
        thetaNewTrans = thetaNew';
        vecThetaNew   = thetaNewTrans(:);
        diff = sum(abs(vecThetaNew-vecTheta))
        vecTheta      = learningRate*vecThetaNew + (1-learningRate)*vecTheta;
    end
    vecThetaOpt    = vecTheta;
else    
    vecThetaOpt    = vecTheta0;
    exitflag       = 0;
end

[~,resEuler] = NLSobjectFunc(vecThetaOpt,setup);
theta   = zeros(setup.nDim1Theta,setup.nDim2Theta);
countIdx = 0;
for i=1:setup.nDim1Theta
    nn = countIdx + sum(setup.select(i,:));
    theta(i,setup.select(i,:) == 1) = vecThetaOpt(countIdx+1:nn)';
    countIdx = nn;
end
nx    = setup.modelPer.nx;
ny    = setup.modelPer.ny;


%% Unfolding theta and accounting for any scaling
% Intercepts - Never scaled
g0             = theta(:,1);    %contains the steady state AND the risk-adj for the controls

% First order terms
scaleX_in_g    = setup.scaleX_in_g;
gx             = theta(:,1+1:1+nx)./repmat(scaleX_in_g',ny,1);

% Nonlinear terms
gxx           = zeros(ny,nx,nx);
gxxx          = zeros(ny,nx,nx,nx);
idx           = 0;
if setup.appOrder > 1
    for i=1:nx
        for j=i:nx
            idx = idx + 1;
            gxx(:,i,j) = theta(:,1+nx+idx)/(scaleX_in_g(i,1)*scaleX_in_g(j,1));
            gxx(:,j,i) = gxx(:,i,j);
        end
    end
end
if setup.appOrder > 2
    idx   = 0;
    for alfa1=1:nx
        for alfa2=alfa1:nx
            for alfa3=alfa2:nx
                idx = idx + 1;
                gxxx(:,alfa1,alfa2,alfa3) = theta(:,1+nx+nx*(nx+1)/2+idx)/(scaleX_in_g(alfa1,1)*scaleX_in_g(alfa2,1)*scaleX_in_g(alfa3,1));
                
                % Using symmetry for alfa1 and alfa2
                if alfa1 == alfa2 && alfa2 ~= alfa3 %alfa1==alfa2~=alfa3
                    %gxxx(:,alfa1,alfa1,alfa3)= gxxx(:,alfa1,alfa2,alfa3);
                    gxxx(:,alfa1,alfa3,alfa1) = gxxx(:,alfa1,alfa2,alfa3);
                    gxxx(:,alfa3,alfa1,alfa1) = gxxx(:,alfa1,alfa2,alfa3);                
                end
                % Using symmetry for alfa2 and alfa3            
                if alfa1 ~= alfa2 && alfa2 == alfa3  %alfa1~=alfa2==alfa3 
                    %gxxx(:,alfa1,alfa2,alfa2)= gxxx(:,alfa1,alfa2,alfa3);                
                    gxxx(:,alfa2,alfa1,alfa2) = gxxx(:,alfa1,alfa2,alfa3);                
                    gxxx(:,alfa2,alfa2,alfa1) = gxxx(:,alfa1,alfa2,alfa3);     
                end            
                % Using symmetry for alfa1,alfa2, and alfa3            
                if alfa1 ~= alfa2 && alfa1 ~= alfa3 &&  alfa2 ~= alfa3 %alfa1~=alfa2~=alfa3
                    %gxxx(:,alfa1,alfa2,alfa3) = gxxx(:,alfa1,alfa2,alfa3);
                    gxxx(:,alfa1,alfa3,alfa2) = gxxx(:,alfa1,alfa2,alfa3);
                    gxxx(:,alfa3,alfa1,alfa2) = gxxx(:,alfa1,alfa2,alfa3);
                    gxxx(:,alfa3,alfa2,alfa1) = gxxx(:,alfa1,alfa2,alfa3);
                    gxxx(:,alfa2,alfa3,alfa1) = gxxx(:,alfa1,alfa2,alfa3); 
                    gxxx(:,alfa2,alfa1,alfa3) = gxxx(:,alfa1,alfa2,alfa3);                                                
                end                
            end
        end
    end
end

% The states
h0   = zeros(nx,1);
hx   = zeros(nx,nx);
hxx  = zeros(nx,nx,nx);
hxxx = zeros(nx,nx,nx,nx);
mx   = setup.mx;
myx  = setup.myx;
if mx > 0
    % Pure Endogenous states
    error('projectStep.m: cannot handled mx > 0')
end
% recall states are in dev from steady state
if myx > 0
    % Laggend controls as states
    h0(mx+1:mx+myx,1)       = g0(1:myx,1)- setup.gSS(1:myx,1);  %contains risk adjustment for the states 
    hx(mx+1:mx+myx,:)       = gx(1:myx,:);
    hxx(mx+1:mx+myx,:,:)    = gxx(1:myx,:,:);
    hxxx(mx+1:mx+myx,:,:,:) = gxxx(1:myx,:,:,:);
end
% Shocks are linear
hx(mx+myx+1:nx,:) = setup.modelPer.hx(mx+myx+1:nx,:);   

% Summarize projection solution
model.scaleX_in_g = scaleX_in_g;
model.appOrder    = setup.appOrder;
model.params      = setup.params;
model.theta       = theta;
model.vecTheta    = vecThetaOpt;
model.labelx      = setup.labelx;
model.labely      = setup.labely;
model.gSS         = setup.gSS;
model.g0          = g0;
model.gx          = gx;
model.gxx         = gxx;
model.gxxx        = gxxx;
model.Gxx         = reshape(model.gxx,ny,nx^2);
model.Gxxx        = reshape(model.gxxx,ny,nx^3);
model.hSS         = setup.hSS;
model.h0          = h0;
model.hx          = hx;
model.hxx         = hxx;
model.hxxx        = hxxx;
model.Hxx         = reshape(model.hxx,nx,nx^2);
model.Hxxx        = reshape(model.hxxx,nx,nx^3);
model.ny          = ny;
model.nx          = nx;
model.nx1         = setup.nx1;
model.myx         = setup.myx;
model.mx          = setup.mx;
model.ne          = size(setup.eta,2);
model.eta         = setup.eta;
model.heteroShocks= setup.modelPer.heteroShocks;
if setup.endoMeanGrid == 1
    varX     = reshape((eye(nx^2)-kron(hx,hx))\reshape((setup.eta*setup.eta'),nx^2,1),nx,nx);
    meanX    = (eye(nx)-hx)\(h0+reshape(hxx,nx,nx^2)*varX(:));
    model.xGrid       = setup.xGrid+repmat(meanX',size(setup.xGrid,1),1);
else
    model.xGrid       = setup.xGrid;
end
model.resEuler    = resEuler;
model.Q           = resEuler(:)'*resEuler(:);
model.exitflag    = exitflag;
model.MAE         = mean(abs(resEuler(:)));
model.RMSE        = sqrt(mean(resEuler(:).^2));
end

